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We discuss an algorithm for the numerical evaluation of NLO multiparton processes. We focus 
hereby on the virtual part of the NLO calculation, i.e. on evaluating the one-loop integration nu- 
merically. We employ and extend the ideas of the subtraction method to the virtual part and we 
use subtraction terms for the soft, collinear and ultraviolet regions, which allows us to evaluate 
the loop integral numerically in four dimensions. A second ingredient is a method to deform the 
integration contour of the loop integration into the complex plane. The algorithm is derived on 
the level of the primitive amplitudes, where we utilise recursive relations to generate the corre- 
sponding one-loop off-shell currents. We discuss the numerical behavior of the approach and the 
application to the leading colour contribution in e + e~ —> n jets, with n up to seven. 
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1. Introduction 

Calculating efficiently multiparton QCD amplitudes and collider observables at next-to-leading 
order (NLO) accuracy is a rather involved task. The most challenging piece is the virtual part. We 
discuss here a numerical approach for the calculation of the virtual part, where we employ the sub- 
traction method and contour deformation [jl|-|n|]. In this regard our algorithm is different from the 
commonly used approaches, based on cut techniques and generalised unitarity or on more tradi- 



tional Feynman graph approaches []12|-|23|], and shows promising features for the implementation 
in a numerical program. The algorithm consists of local subtraction terms to subtract divergences 
arising from the soft, collinear and ultraviolet (UV) regions of the virtual part, which render the 
integrand finite in the respective regions, and of a method to deform the integration contour of the 
loop integration into the complex plane in order to circumvent the remaining on-shell singularities. 
It works on the level of colour-ordered primitive amplitudes, where we utilise recursive algorithms 
to compute the corresponding one-loop off-shell currents for the bare primitive amplitudes, and 
is therefore fast and easily implemented. The numerical loop integration is performed together 
with the integration over the phase-space of the external particles in one Monte Carlo integration. 
The subtraction terms yield simple results upon analytic integration over the loop-momentum, and 
the resulting pole structures cancel exactly against the pole structures from the soft and collinear 
parts of the real emission contributions as well as of the UV counterterm from renormalisation. 
The algorithm goes hand in hand with the usual subtraction method for the real emission contri- 
butions, where we employ Catani-Seymour dipole subtraction [|4]-|^]. We applied the method to 
e + e~ — > n jets, with n up to seven, in the large-Af c limit. Up to four jets we reproduce the known 
results for the respective jet rates with very good agreement. Increasing the number of jets up to 7 
shows a good scaling behaviour in CPU time with respect to the number of final state particles. 

2. Subtraction method 

The contributions to an infrared-safe observable at next-to-leading order with n final state 
particles can be written in a condensed notation as 

(Of LO = J O n+l da R + J O n do v + J O n da c , (2.1) 

n+\ n n 

where do R denotes the real emission contribution, which corresponds to the square of the Born am- 
plitude with (n + 3) partons |=2^v 3 | 2 , da v denotes the virtual contribution, which corresponds to the 
interference term of the renormalised one-loop amplitude with the Born amplitude 2D\e(£/j+2 ^+2)' 
and da c subtracts initial state collinear singularities. Each term is separately divergent and only 
their sum is finite. One adds and subtracts suitably chosen pieces to be able to perform the phase- 
space integrations and the loop integration by Monte Carlo methods: 

{0) nlo = J ( 0n+l do R -O n do A ) + J {O n da^ e -O n da L ) 

n+\ n+loop 
+ J (OndO^j + O n J dO L + O n j dO A + O n dO C ^j , (2.2) 



loop 
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where we write the renormalised virtual piece as the sum of the bare part and a counterterm part 

J O n da v = Jo n J dal ie + J O n do v CT . (2.3) 

n n loop n 

The first term (O ll+ ido R — O n do A ^ is by construction integrable over the (n+ l)-particle phase- 
space and can be evaluated numerically. In the dipole subtraction the subtraction term O n do A is 
given by a sum over dipoles. The second term (O n dc^ we — O n do L ) is by construction integrable 
over the rc-particle phase-space and the loop-momentum space in four dimensions. We thus extend 
the subtraction method to the virtual part such that we can evaluate the integral of the one-loop 
amplitude numerically. In any Monte Carlo integration the error scales independently of the di- 
mensionality of the integration region. The phase-space and loop integrals can thus be evaluated in 
a combined Monte Carlo integration and no additional costs to evaluate the loop integral separately 
per phase-space point are introduced. After analytical integration of the subtraction terms over the 
unresolved one-parton phase-space and the loop-momentum space respectively, the IR poles of the 
virtual subtraction term cancel with the IR poles of the real subtraction term and the initial state 
collinear subtraction term, whereas the UV poles of the virtual subtraction term cancel with the UV 
poles of the counterterm. Therefore the third term is also finite and can be evaluated numerically. 
In short, the NLO contribution is given as the sum of three finite contributions 

(Of LO = (O)fif + (O)^ + <« tlon . (2.4) 



3. Colour decomposition and kinematical setup 

Amplitudes in QCD may be decomposed into group-theoretical colour factors multiplied by 
purely kinematical partial amplitudes: 

^ (1) =£CA (1) (3-1) 

i 

In the colour-flow basis [28, 3^, the colour structures are linear combinations of monomials in 
Kronecker 5,/s. One-loop partial amplitudes are further expressed in terms of linear combinations 
of primitive one-loop amplitudes, where a primitive amplitude is defined as a gauge-invariant set of 
colour-stripped Feynman diagrams with a fixed cyclic ordering of the external partons and a definite 



routing of the external fermion lines through the diagrams [32], which are important properties for 
our method. We will drop any subscripts refering to colour on the primitive amplitudes A' 1 ) from 
now on. Due to the fixed cycling ordering there are only n different loop-propagators occuring in 
a primitive amplitude with n external legs. With the notation as in fig.([j]) we define kj = k — qj, 
with qj = YJi=\ Pu where k denotes the loop-momentum in clock-wise direction, and the pi are the 
outgoing external momenta. We can write the bare primitive one-loop amplitude, in dimensional 
regularisation, as 

^=/(0Xi> effi. = iWJJ^z^a. (3 - 2) 

where P a (k) is a polynomial of degree a in the loop momenta k and the +/5-prescription tells 
us in which direction the poles of the propagators should be avoided. Soft singularities arise for 
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Pn-1 




Figure 1: The labelling of the momenta for a primitive one-loop amplitude. 

k ~ qi and pj = mj_ l , mj = 0, pj +1 = mj +1 , i.e. a massless particle exchanged between two on-shell 
particles and soft k { . Collinear singularities arise for k ~ q t —xpi and pj = 0, m^_i = 0, m,- = 0, 
where x £ [0, 1], i.e. a massless external on-shell particle attached to two massless propagators and 
collinear momenta ki-\ , kj and UV singularities arise for \k\ — > oo and 4 + a — 2n > 0. 

4. Local infrared and ultraviolet subtraction terms 

The soft and collinear subtraction terms for massless QCD read 

fl.W -AjY Pj-Pj+l c W(lJ2 1.2 7,2 nn.(0) 

w t K j-i K j K j+i 



7,2 7,2 ' 7,2 7,2 

j— i j K j K i+i 



Af\ (4.2) 



where the sum over j € I g is over all gluon propagators j inside the loop. Furthermore, Sj = 1 if 
the external line j corresponds to a quark and Sj = 1/2 if it corresponds to a gluon. The functions 
<?soft an d £con ensure tnal: the integration over the loop momentum is UV finite. They are chosen 
to be equal to one in the corresponding soft and collinear limits and to suppress the integrands 
by appropriate powers of l/\k\ in the ultraviolet limit \k\ — > °°. The soft and collinear subtraction 
terms yield simple results upon analytic one-loop integration in dimensional regularisation, with 
the appropriate pole structure to cancel the poles from the real emission and initial state collinear 
subtraction. We derived them for the massless as well as the massive case. 

The UV subtraction terms correspond to local counterterms g^ jW (k,Q,{pj},{mj}) to propa- 
gator and vertex corrections G. These local counterterms are used in a recursive approach to build 
a total UV subtraction term to the bare primitive amplitude. The local counterterms are obtained 
by expanding the relevant loop propagators around a new UV propagator (k 2 — /i^y) _1 , where 
k = k — Q and for a single propagator we have 

i _ i . + 21 ■ (p - Q) _ (g - Of + ^ + [a ■ (p - 0)] 2 + e IX y 



(k-pf P-^v (P-figv) 2 (P-«5v) 2 (P-^v) 3 VI*! 5 

We can always add finite terms to the subtraction terms. For the local counterterms we choose 
the finite terms such that the finite parts of the integrated local counterterms are independent of Q 



4 



Numerical evaluation ofNLO multiparton processes 



C. Reuschle 



and proportional to the pole part, with the same constant of proportionality. The integrated local 
counterterms assume then the general form 

C G ' (0) ({^},H})(i-ln^) + ^( £ ) (4.4) 

where C G ^°'({pj},{nij}) takes on the typical form of the corresponding renormalised colour- 
ordered Born level propagator or vertex function. This ensures that the sum of all integrated local 
counterterms is again proportional to a tree-level amplitude. We derived the local UV counterterms 
for the massless as well as the massive case. 



5. Contour deformation 

Having a complete list of local UV and IR subtraction terms at hand, we can ensure that the 
integration over the loop-momentum gives a finite result and can therefore be performed in four 
dimensions. However, there is still the possibility that some of the loop-propagators go on-shell 
for real values of the loop-momentum. Therefore we shift the integration contour into the complex 
space C 4 , where the integration contour must be chosen such that whenever possible the poles of 
the propagators are avoided. We set 

k = k + iK(k), (5.1) 
where contains only real components. After the deformation our one -loop integral reads 



d A k R(k) f d A k 



R(k(k)) 

(5.2) 



n [kj-mj-K 2 + 2ikj-K 



7=1 v J 7=1 

where we integrate over the four real components in h 1 . To match Feynman's +z5-prescription 
we have to construct the deformation vector K such that k; ■ fc > whenever k f — m 2 - = 0, and the 

* J J 

equal sign applies only if the contour is pinched. If the contour is pinched the singularity is either 
integrable by itself or there is a subtraction term for it. The numerator function R{k) has only poles 
at k — /luv = 0- Choosing jU^v sufficiently large on the negative imaginary axis ensures that the 
integration contour always stays away from these poles. The form of the deformation vector k has 
a direct impact on the Monte Carlo integration error. We employ several techniques to reduce the 
Monte Carlo integration error. First, we split our integral into an exterior and an interior part: 

r d 4 k r d 4 k r d 4 k " k 2 --m 2 

This splitting is holomorphic in k, which means we can use different contours for the evaluation 
of the two parts. In the exterior part I ext the poles from k 2 — m 2 are absent and the deformation 
vector K can simply be chosen as = g^ v (k v — <2 V ) • If we choose the arbitrary four- vector Q in 
k = k — Q to be 

n 7=1 
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we can arrange that the integrand of the interior part drops off with two extra powers in \/\k\ 
for \k\ — > oo and thus receives less contributions from the UV region. The construction of the 
deformation vector for the internal region is along the lines of [Q] and given in detail in 

The lines connecting the origins of the light cones given by (k — qj) = in fig.@ are the 
regions of collinear singularities. If the loop-momentum approaches a collinear singularity we have 
a subtraction term for it. To improve the numerical behaviour in the vicinity of such a collinear 
singularity we further split the interior part into several sub-channels, where each sub-channel 
corresponds to one line segment in fig.(EJ). We write 

l^ m= tJik w ' (k)m <5 - 5) 

n-l 

with w,- > and £ w, (£) = 1. For the weights w,- we use 

i=\ 

7=1 V 7 

with a = —2. This division into sub-channels is not holomorphic and we have to use the same 
integration contour in all the sub-channels. However, we can use different parametrisations, i.e. 
coordinate systems, for the different sub-channels. A possible choice to parametrise a certain sub- 
channel is the generalisation of elliptical coordinates to four dimensions. The coordinate system 
of the sub-channel j is thereby chosen such that the origin of this coordinate system corresponds 
to the line segment, which connects qj with qj+i. This is shown in fig.(||). The sampling for an 
individual sub-channel is discussed in detail in [0J. Upon this parametrisation we observe certain 
periodic oscillations of the integrand in the generalised elliptical coordinates. In order to average 
out these oscillations we use the method of antithetic variates, as described in [01, which improves 
the numerical behaviour significantly. 

6. Numerical stability in the UV regions 

It turns out that also the UV regions give rise to large numerical oscillations. In their stan- 
dard form the local UV counterterms are constructed such that after subtraction we are left with a 
leading l/|&| 5 -behaviour for \k\ — > °°, which is formally enough to ensure UV fmiteness. In order 
to improve the damping in the UV region, we choose the expansion of the UV denominators such 
that we are left with a leading l/|£| 7 -behaviour, where our experience shows that we only need to 
do this for the corrections to propagators and three-valent vertices. We do the same for the UV- 
suppression functions g^ t and in the soft and collinear subtraction terms. Together with the 



choice for Q in eq.(5.4) the UV-behaviour of the integrand is significantly enhanced. More details 



and additional comments can be found in [Eh. 

7. Recursion relations 

We use Berends-Giele type recursion relations [ |33| ] to compute the tree amplitude, the bare 
one-loop integrand G^ IS and the total UV subtraction term Guy. These recursion relations are 
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Figure 2: Diagram (a) shows the origins of the light 
cones in electron-positron annihilation. The origins are 
given by the (n — 1) vertices. The vertices are con- 
nected by (n — 2) line segments. We decompose /,„, into 
[n — 2) sub-channels, such that each sub-channel corre- 
sponds to one line segment. This is shown in diagram 
(b), where the line segment from q\ to q2 is emphasised. 



Figure 3: The diagram shows the coordinate sys- 
tem for the y'th sub-channel. We use generalisations 
of elliptical coordinates to four dimensions. At the 
origin the radial component of the loop-momentum 
in the elliptical parametrisation equals zero and the 
loop -momentum coincides with a point on the line 
segment that connects qj with qj+i- 




n- 1 




+ 





Figure 4: Recursive relations for a three-valent toy model. 



shown in fig.@ for the case of a three-valent toy model. The recursive relations have been im- 
plemented for all necessary cases, i.e. gluon currents as well as quark and antiquark currents, that 
are needed for e + e~ — > n jets in the large-./V c limit. In the direct loop contributions, where the off- 
shell leg couples directly through a vertex to the loop, the two edges of the vertex are connected 
to loop-propagators. We can cut open one of these two edges by replacing the tensor structure 
of the corresponding propagator by a sum over (pseudo)-polarisations. More details can be found 
in [0| . We can check the cancellations in the UV region between the recursive constructions of 
the bare one-loop integrand and the total UV subtraction term. This is shown in fig.(|5]), where 
we plot \2<Rz(A(°> G' ))| vs. a UV scaling parameter Auv for e + e — > 3 jets, in leading colour 
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ie+06 




1 10 100 



Figure 5: The plot shows the helicity summed \2U\c(A^*G^)\ vs. a UV scaling parameter Auv for 
e + e~ — > 3 jets, in leading colour approximation, where we scale a fixed value for the loop-momentum ac- 
cording to k = k — <2fi X ed = Auv^fixed- The unsubtracted total integrand (upper fit in grey) falls off locally like 
1/|£| 3 , which leads to UV divergences upon integration. The (standard) UV subtracted total integrand (mid- 
dle fit in blue) falls off like 1/|&| 5 , which is clearly UV finite. The (numerically improved) UV subtracted 
total integrand (lower fit in pink), which contains those local UV counterterms that also subtract the 1/|£| 5 - 
and 1/ |^| 6 -behaviour, falls off like l/\k\ 7 . 



approximation. 



8. NLO results for n jets in electron-positron annihilation 



We have calculated the Durham jet rates in electron-positron annihilation in the leading colour 
approximation for up to seven jets The cross section for n jets normalised to the leading order 
(LO) cross section for e + e~ — > hadrons reads 



n-2 



a„(m) + 



a,(ju) 



n-\ 



(Join) \ 2k J " n ^ r ~> ' y 2% 
One can expand the perturbative coefficient A n and B n in l/N c : 

n-2 



B n (n) + 0(cc?). 



(8.1) 



A n =N C 



N c 



B n = N c 



n-\ 



B n i c + 



(8.2) 



We calculate the LO coefficient A„ / c and the NLO coefficient B n j c for n < 7 at the renormalisation 
scale jU equal to the centre of mass energy. We take the centre of mass energy to be equal to the 
mass of the Z-boson. The scale variation can be restored from the renormalisation group equation. 
The calculation is done with five massless quark flavours. Fig.@ shows the comparison of our 
numerical approach with the well-known results for two, three and four jets [[34], 35]. We observe 
an excellent agreement. The results for five, six and seven jets for a jet parameter of y cut = 0.0006 
are shown in fig.(0). Fig.(||) shows the CPU time required for one evaluation of the Born contri- 
bution, the insertion term and the virtual term as a function of the number n of jets. One notes 
that the insertion term is almost as fast as the Born contribution. For all contributions the CPU 
time per evaluation increases only very moderately as a function of n. Within our method all three 
contributions scale asymptotically as n 4 , as expected [36]. The virtual part has the same scaling 



behaviour as the Born contribution. This moderate growth imposes almost no restrictions on the 
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Vm ym yaa 



Figure 6: Comparison of the NLO corrections to the two-, three- and four-jet rate between the numerical 
calculation and an analytic calculation. The error bars from the MC integration are shown and are almost 
invisible. 

CPU time 

T| N c {%) n - 2 A nJ c I A^frX/c 

5 (2.4764 ± 0.0002) • 10 4 (1.84±0.15) • 10 6 

6 (2.874±0.002)-10 5 (3.88 ±0.18) • 10 7 

7 (2.49 ± 0.08) • 10 6 (5.4 ± 0.3) • 10 8 

Figure 7: The results for the LO and NLO leading- 
colour jet rate coefficients, for a jet parameter of y cut = Fi gure 8: CPU time required for one evaluation of 
0.0006. Shown are the results for five, six and seven the various contributions as a function of the number 
j e t s n of jets. The times are taken on a single core of a 

standard PC. 

number of final state partons to which our method can be applied. The practical limitations arise 
from the fact that the number of evaluations required to reach a certain accuracy increases with n. 
This behaviour is already present at the Born level and not inherent to our method. The calculation 
of the seven-jet rate takes about five days on a cluster with 200 cores. 
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